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SUMMARY 

The  spectral  element  method  for  the  two-dimensional  shallow  water  equations  on  the  sphere  is  presented. 
The  equations  are  written  in  conservation  form  and  the  domains  are  discretized  using  quadrilateral 
elements  obtained  from  the  generalized  icosahedral  grid  introduced  previously  (Giraldo  FX.  Lagrange  - 
Galerkin  methods  on  spherical  geodesic  grids:  the  shallow  water  equations.  Journal  of  Computational 
Physics  2000;  160:  336-368).  The  equations  are  written  in  Cartesian  co-ordinates  that  introduce  an 
additional  momentum  equation,  but  the  pole  singularities  disappear.  This  paper  represents  a  departure 
from  previously  published  work  on  solving  the  shallow  water  equations  on  the  sphere  in  that  the 
equations  are  all  written,  discretized,  and  solved  in  three-dimensional  Cartesian  space.  Because  the 
equations  are  written  in  a  three-dimensional  Cartesian  co-ordinate  system,  the  algorithm  simplifies  into 
the  integration  of  surface  elements  on  the  sphere  from  the  fully  three-dimensional  equations.  A  mapping 
(Song  Ch,  Wolf  JP.  The  scaled  boundary  finite  element  method — alias  consistent  infinitesimal  finite 
element  cell  method — for  diffusion.  International  Journal  for  Numerical  Methods  in  Engineering  1999;  45: 
1403-1431 )  which  simplifies  these  computations  is  described  and  is  shown  to  contain  the  Eulerian  version 
of  the  method  introduced  previously  by  Giraldo  ( Journal  of  Computational  Physics  2000;  160:  336-368) 
for  the  special  case  of  triangular  elements.  The  significance  of  this  mapping  is  that  although  the  equations 
are  written  in  Cartesian  co-ordinates,  the  mapping  takes  into  account  the  curvature  of  the  high-order 
spectral  elements,  thereby  allowing  the  elements  to  lie  entirely  on  the  surface  of  the  sphere.  In  addition, 
using  this  mapping  simplifies  all  of  the  three-dimensional  spectral-type  finite  element  surface  integrals 
because  any  of  the  typical  two-dimensional  planar  finite  element  or  spectral  element  basis  functions 
found  in  any  textbook  (for  example,  Huebner  et  al.  The  Finite  Element  Method  for  Engineers.  Wiley,  New 
York,  1995;  Karniadakis  GE,  Sherwin  SJ.  Spectraljhp  Element  Methods  for  CFD.  Oxford  University 
Press,  New  York,  1999;  and  Szabo  B,  Babuska  I.  Finite  Element  Analysis.  Wiley,  New  York,  1991)  can 
be  used.  Results  for  six  test  cases  are  presented  to  confirm  the  accuracy  and  stability  of  the  new  method. 
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1.  INTRODUCTION 

The  spectral  element  method  combines  the  benefits  of  both  the  spectral  and  finite  element 
methods.  The  spectral  element  method  can  automatically  generate  any  order  polynomial  basis 
function,  as  in  the  spectral  method,  while  allowing  for  the  geometrical  flexibility  enjoyed  by  the 
finite  element  method.  This  geometric  flexibility  permits  the  use  of  an  unstructured  grid.  The 
spectral  element  method  may  use  any  type  of  Jacobi  polynomial  to  define  the  basis  functions, 
but  typically  either  Chebyshev  or  Legendre  polynomials  are  used;  in  this  paper,  Legendre 
polynomials  are  used.  The  polynomial  order  of  the  spectral  element  method  can  be  generated 
automatically  along  with  the  corresponding  numerical  integration  rule,  which  makes  the 
method  very  suitable  for  constructing  general  codes.  By  selecting  the  numerical  integration 
(quadrature)  points  to  be  the  collocation  points  results  in  a  diagonal  mass  matrix  that  reduces 
the  storage  requirements.  The  efficiency  of  the  method  is  further  enhanced  by  not  having  to 
define  the  basis  functions  explicitly  because  implicit  relations  for  the  inner  products  of  the 
functions  and  their  derivatives  can  be  defined  a  priori  [1].  Since  the  collocation  points  are  not 
equi-spaced,  staggered  grids  can  be  generated  automatically  by  using  varying-order  polynomi¬ 
als  for  the  different  variables  (say  the  pressure  and  velocity  in  Navier-Stokes)  as  is  done  in 
Reference  [2]  for  the  shallow  water  equations.  This  satisfies  the  Babuska-Brezzi  (or  inf-sup) 
condition,  thereby  avoiding  the  development  of  spurious  pressure  modes.  Many  researchers 
have  used  the  spectral  element  method  successfully  for  advection- diffusion  [3],  Stokes  flow  [1], 
Navier-Stokes  [4],  and  the  shallow  water  equations  [2,5,6]. 

The  recent  paradigm  shift  in  large-scale  computing  from  vector  machines  (having  few  but 
very  powerful  processors)  to  distributed  memory  machines  (having  many  but  less  powerful 
processors)  has  led  atmospheric  scientists  to  explore  methods  other  than  the  spectral  method 
for  solving  the  governing  equations  on  the  sphere.  (The  shallow  water  equations  have  been 
used  customarily  to  develop  new  numerical  methods  for  atmospheric  models  because  they 
exhibit  the  same  wave  behavior  as  the  more  complex  baroclinic  equations  governing  the 
motion  of  the  atmosphere.)  Besides  not  parallelizing  well,  the  spectral  method  also  suffers  from 
the  restriction  that  the  grid  be  a  longitude-latitude  grid,  which  packs  too  many  unnecessary 
points  at  the  poles.  By  exploring  other  classes  of  methods,  researchers  are  also  free  to  choose 
other  types  of  grids.  For  example,  in  Reference  [6]  a  cubic  gnomonic  grid  is  used  in 
conjunction  with  the  spectral  element  method.  In  Reference  [7]  an  icosahedral  grid  is  employed 
with  the  finite  difference  method.  In  Reference  [8]  a  spiral  triangular  grid  similar  to  the 
icosahedral  grid  is  used  with  the  weak  Lagrange-Galerkin  finite  element  method.  (For  a 
complete  review  of  grids  for  tiling  the  sphere  see  Reference  [9].)  However,  many  of  these  new 
approaches  continue  to  follow  in  the  footsteps  of  the  spectral  method  by  writing  the  equations 
in  spherical  co-ordinates;  the  only  exception  is  the  finite  difference  method  presented  in 
Reference  [7],  which  solves  the  equations  in  a  co-ordinate  invariant  form. 

Although  spherical  co-ordinates  seem  to  be  the  natural  choice  on  the  sphere,  they  present 
many  problems  and  associated  computational  costs.  As  mentioned  above,  the  spectral  method 
on  the  sphere  requires  the  use  of  a  longitude -latitude  Gaussian  grid,  which  introduces  too 
many  redundant  points  around  the  poles  and  this  situation  is  exacerbated  as  the  resolution 
increases.  Finite  difference,  finite  element,  and  spectral  element  methods  suffer  from  the 
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existence  of  singularities  in  the  equations  at  the  poles  in  this  co-ordinate  system.  This  difficulty 
can  be  circumvented  by  applying  rotation  transformations  (as  is  done  in  Reference  [6]  from 
spherical  to  gnomonic  space)  or  by  using  Cartesian  co-ordinates. 

In  References  [10-13]  the  equations  are  solved  in  Cartesian  rather  than  spherical  co¬ 
ordinates.  In  Reference  [10]  the  Lagrange  multiplier  approach  for  transforming  the  shallow 
water  equations  from  spherical  to  Cartesian  co-ordinates  is  introduced.  This  idea  is  then  used 
to  construct  a  semi-Lagrangian  shallow  water  model  but  on  a  longitude -latitude  grid,  an 
unnecessary  remnant  of  spectral  methods.  In  Reference  [13]  the  Taylor -Galerkin  method  in 
three-dimensional  Cartesian  space  using  an  icosahedral  grid  is  presented.  However,  in  Refer¬ 
ence  [13],  the  work  stops  short  by  then  mapping  the  three-dimensional  linear  triangles  to  a 
local  two-dimenaional  space;  an  unnecessary  step  since  all  the  computations  can  be  carried  out 
in  three-dimensional  space,  thereby  avoiding  any  type  of  mapping  to  a  local  space.  In 
References  [11,12]  we  developed  a  weak  Lagrange -Galerkin  finite  element  method  in  three- 
dimensional  Cartesian  space  using  triangular  linear  basis  functions  on  an  icosahedral  grid;  this 
method  can  only  use  linear  elements,  however. 

In  the  current  paper,  the  equations  are  once  again  formulated  in  three-dimensional  Cartesian 
space  (as  in  References  [11,12]),  thus  simplifying  the  construction  of  surface  integrals  directly 
onto  the  three-dimensional  space.  One  of  the  many  differences  between  this  work  and  that 
presented  in  References  [11,12]  is  that  it  permits  the  use  of  very  high-order  quadrilateral 
spectral  elements  (up  to  32nd  order)  rather  than  linear  triangular  elements.  Quadrilaterals  had 
to  be  adopted  over  triangles  because  no  C°  nodal  expansion  exists  for  triangles  (modal 
expansions  do  exist  however,  see  References  [14,15]).  Some  promising  ideas  for  constructing 
spectral  elements  on  triangles  can  be  found  in  References  [16,17]  and  we  hope  to  return  to  this 
in  a  future  paper. 

In  order  to  be  able  to  construct  these  high-order  elements  directly  onto  the  three- 
dimensional  Cartesian  space,  a  special  mapping  had  to  be  developed  (this  mapping  was  derived 
independently  by  the  current  author  and,  long  before,  by  Song  and  Wolf  [18]).  In  order  to  test 
whether  the  mapping  is  mathematically  correct,  the  Eulerian  method  described  in  References 
[11,12]  was  compared  against  the  new  mapping  for  the  special  case  of  linear  triangular 
elements;  in  Appendix  A  we  show  that,  for  this  case,  the  two  methods  yield  equivalent  results 
thereby  proving  that  the  mapping  is  indeed  correct. 

Writing  the  equations  in  Cartesian  space  and  using  the  mapping  described  herein  allows  the 
discretization  of  the  equations  in  conservation  form;  this  has  beneficial  implications  for  the 
conservation  properties  of  the  scheme  and  it  also  yields  a  form  much  more  suitable  for  using 
standard  computational  fluid  dynamic  (CFD)  techniques.  In  fact,  in  this  form  the  shallow 
water  equations  appear  as  a  non-homogeneous  case  of  the  Euler  equations  that  govern  the 
motion  of  inviscid  compressible  flow. 

The  general  triangular  icosahedral  grid  introduced  in  Reference  [12]  is  used  for  constructing 
the  higher-order  spectral  elements  on  the  sphere.  Each  triangle  is  divided  into  three  quadrilat¬ 
erals  because  we  require  quadrilateral  elements  in  order  to  construct  a  high-order  C°  nodal 
expansion.  The  current  paper  represents  a  departure  from  previously  published  work  on 
solving  the  shallow  water  equations  on  the  sphere  in  that  the  equations  are  all  written, 
discretized,  and  solved  in  three-dimensional  Cartesian  space,  which  simplifies  all  of  the 
calculations.  This  simplification  is  especially  useful  when  dealing  with  high-order  spectral 
element  basis  functions  on  a  curved  surface  (i.e.,  a  sphere). 
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The  method  described  herein  can  be  used  to  construct  any  type  of  high-order  finite  element 
method,  including  the  nodal  expansion  described  here  and  the  modal  expansion  using 
triangular  or  quadrilateral  hierarchic  basis  functions  described  in  References  [14,19].  However, 
this  paper  only  focuses  on  a  spectral-type  finite  element  method  (nodal  expansion)  on 
quadrilateral  grids.  In  the  next  section,  we  describe  the  governing  equations  used  in  this  paper, 
namely,  the  shallow  water  equations. 


2.  SHALLOW  WATER  EQUATIONS 


The  two-dimensional  spherical  shallow  water  equations  in  Cartesian  conservation  form  are 
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but  note  that  if  we  move  the  pressure  terms  to  the  right-hand  side,  we  get 
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This  system  can  now  be  written  more  compactly  as 
d(p 

^  +  V-(pu)  =  S  (?) 
where 
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where  Q  and  a  are  the  rotation  of  the  earth  and  its  radius  respectively.  In  the  above  relations 
FP  is  the  force  due  to  pressure,  Fp  is  the  force  due  to  the  rotation  of  the  earth  (Coriolis),  and 
Fc  is  the  force  required  to  constrain  the  fluid  particles  to  remain  on  the  surface  of  the  sphere, 
and  /<  is  a  Lagrange  multiplier  used  to  satisfy  this  condition.  This  term  is  best  obtained  in 
discrete  form.  If  it  is  obtained  from  the  continuous  equations,  then  for  the  discrete  equations 
the  particles  may  no  longer  be  guaranteed  to  remain  on  the  surface  of  the  sphere. 

The  equations  are  solved  for  the  four  conservation  variables  q>  =  (cp,  cpu,  <pv,  q>w),  where  cp 
denotes  the  geopotential  height  and  u  the  velocity  field.  Clearly,  we  could  also  include  other 
forcing  functions  and  they  would  then  be  included  in  S (<p).  In  the  next  section,  we  show  how 
to  discretize  Equation  (2). 


3.  DISCRETIZATION 


3.1.  Spatial  discretization 

Beginning  with  Equation  (2)  we  can  define  a  spectral  element  formulation  by  taking  the  weak 
form 


Jq 
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dQ  =  0 


and  integrating  by  parts  (Green’s  theorem)  such  as 
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where  it  is  understood  that  the  integrals  with  subscripts  Q  denote  area  integrals  and  those  with 
T  denote  boundary  integrals.  In  the  case  of  flow  on  a  sphere,  the  second  integral  vanishes 
because  the  basis  functions  are  chosen  to  be  continuous  across  element  interfaces  and  there  are 
no  true  boundaries,  in  other  words 


n-ui J/(p  dF  =  0 


and  we  are  left  with 
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Jq  Jq  Jq 

By  letting  the  variables  (p  be  approximated  by  XfT ,  i//,r/)„  where  NP  represents  the  total  number 
of  collocation  points  within  each  element,  we  can  then  write  relations  for  the  integrals  above 
(the  basis  functions  i//  are  described  later  in  this  paper).  This  yields  the  following  integral 
equation: 
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which  can  then  be  written  in  the  matrix  form 
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where  it  is  understood  that  i,j,k=  l, ,  NP  and  so  My  represents  an  NP  x  NP  matrix,  and 
so  on. 
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3.1.1.  Element  matrices.  Substituting  in  the  forces  S (<p)  from  the  right-hand  side  of  Equation 
(1)  into  Equation  (4)  allows  us  to  write  all  of  the  element  matrices  as  follows: 
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The  matrices  M,  A,  /J\  and  1<X  represent  the  mass,  advection,  pressure,  and  rotation  matrices 
respectively.  All  of  these  matrices  are  independent  of  the  conservation  variables  (p  and  only 
need  to  be  computed  once  because  they  remain  constant  throughout  the  integration.  The 
Lagrange  multiplier  /.t  has  not  been  included  in  the  force  vector  S:(<p).  This  term  is  actually 
computed  at  the  discrete  level  and  only  after  each  time  integration  has  been  completed.  The 
procedure  used  to  compute  this  term  is  explained  in  the  section  on  the  time  discretization 
scheme. 

3.1.2.  Basis  functions.  The  conservation  variables  belong  to  the  following  space: 

(tp,  (pu,  <pv,  (pw)eHl(Q) 

and  their  test  functions,  likewise,  are  contained  in  the  space 
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foeHH  Q) 


In  other  words,  they  belong  to  the  set  of  square  integrable  functions  whose  first  derivatives  are 
also  square  integrable  (Sobolev  space).  Note  that  no  boundaries  are  defined  in  this  space  at  all, 
this  is  because  flow  along  the  surface  of  a  sphere  has  no  boundaries.  The  element  basis 
functions  are  written  in  the  following  form  [18]: 


MZ,  >h  0  =  Cte  >/)  (9) 

where  t  denotes  the  outward  pointing  co-ordinate  of  the  surface  defined  by  the  shape  functions 
///,(C,  tj );  t  =  0  gives  the  center  of  the  sphere  and  t  =  1  yields  the  surface  of  the  sphere.  The 
basis  functions  given  by  Equation  (9)  simplify  the  construction  of  high-order  methods  on  the 
sphere  or  any  other  curved  surface  because  the  functions  ^  can  be  the  typical  finite  element  or 
spectral  element  basis  functions  in  R2.  Without  this  type  of  mapping  it  would  not  be  possible 
to  construct  high-order  basis  functions  on  curved  geometries  in  Cartesian  space  as  we  have 
done  in  this  paper. 

The  shape  functions  defining  a  specific  surface  element  are  obtained  as  a  tensor  product 
from  the  one-dimensional  Legendre  cardinal  functions  as 


$,(Z,ti)  =  hj(e)hk(ri)  (10) 

where  i  =  l, ,  N ,,,  j,k=\, . .  .  ,P  and  NP  =  P2.  The  integer  P  denotes  the  number  of 
Legendre-Gauss-Lobatto  (LGL)  points  in  each  direction  (c  and  ;/ ),  and  is  equal  to  P  =p  +  1, 
where  p  denotes  the  polynomial  order  of  the  Legendre  cardinal  functions  (see  References 
[1-3,5]  for  further  details  on  these  basis  functions). 

3.1.3.  Integration.  In  order  to  integrate  the  matrices  defined  by  Equations  ( 5)  (8)  they  must 
first  be  transformed  from  the  global  R3  space  (Q  =  Q(x,  y,  z))  onto  the  the  local  R2  surface 
element  space  (D  =  Q(c, ;/)). 

The  Jacobian  of  the  transformation  is  given  by 


dx 

dx 

dx 

dt] 

~dC 

dy 

dy 

dy 

dr/ 

dt 

dz 

dz 

dz 

dt] 

dt 

with  the  inverse  Jacobian  being 
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dd_ 

dd 

dd 

dx 

dy 

dz 

dr] 

dr ] 

dr] 

dx 

8y 

dz 

dd 

dd 

dd 

dx 

dy 

dz 

which  in  terms  of  the  physical  co-ordinates  is 


J-' 


1  /  dy  dr\  l  (  dz 

d\dr]"  y  dr])  1  [)  dd 
if  dz  fe  \  if dx 

~c{xdf,~8f,z) 

1  f dx  djA  if  dy 

t{df,y-xdf,) 


dy  A  l  dy  dz  dy 
dd)  \dd  di]  di]dd) 
dz  \  /  dx  dz  dx  dz  \ 

X  dd  J  \drj  dd  dd  drj  J 

dx  \  /  dx  dy  dx  dy\ 

dd /  \d£  dr\  dr]  dd  J 


(12) 


The  determinant  of  the  Jacobian  matrix  (11)  can  be  written  in  the  following  form: 


\J  I 


dx 


•G 


where 


dx  dx 


(13) 


and  |G|  denotes  the  magnitude  of  the  normal  vector  to  the  surface  element  Q  in  local  space 
[20].  All  of  these  terms  are  easily  obtained  from  the  basis  functions  t//  and  the  spatial 
co-ordinates  x  from  the  approximations 

Np 

x=  X 

i  =  1 

Having  defined  the  mapping  from  physical  to  computational  space,  we  can  now  integrate 
the  element  matrices  numerically.  The  mass,  advection,  pressure,  and  rotation  matrices  become 


Mjj  = 


/fx  — 
Ai]k  — 


\ I'i'I'j  dH  =  ^  wqd2\G(dq,  t]q)\\j/i(dq,  tlgWjitq,  >1q) 


<7  =  1 


dfx  Wk  d-Q=  I  wqd2\G(dq,  >lq)\( d^^+d  ?§! ’  +  &  }u) 


q=i 


dd  dx  dr]  dx 


dx 
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px  _ 
'  ijk  — 


,  /  #/fcHO  ^  MTr;  up  J  r^k^  r8ij/k8tj  f  8C 

M  dQ  =  X  w’/  lGU?>  >/,)|^(^  m f  tit  ft 


l 

9=  1 


dx  dt]  dx  8x 


R?jk  = 


'l'i'1'j'l'k'l'l-fk  i  WqC4\ |G(^,  tjq)\lj}i(^q,  riq)\j/j(Zq,  l1q)']'k(Zq,  *1  q)Ml(£q,  flq) 

9=1 


X, 

Cl 


where 


A  ijk  =  A^.i  +  A^j  +  A^.k 

and  i,  j,  k  denote  the  Cartesian  directional  vectors.  The  value  of  C  need  not  be  of  any  concern 
because  at  the  surface  of  the  sphere  it  is  defined  by  C  =  1  •  The  integer  Ne  represents  the  total 
number  of  LGL  quadrature  points  required  within  each  element. 

The  number  of  quadrature  points  used  are  NQ  =  Q2,  where  Q  represents  the  number  of  LGL 
quadrature  points  in  each  direction  (c  and  if).  In  general,  for  exact  integration  of  the  matrices 
we  require  Q  =  (cP  +  l)/2  quadrature  points,  where  c  is  an  integer  constant  denoting  the  factor 
of  the  maximum  order  matrix.  For  the  spherical  shallow  water  equations  c  =  4,  corresponding 
to  the  rotation  matrices  Rx,  thereby  requiring  Q  =  (4 P  +  l)/2  quadrature  points.  Numerically, 
it  was  found  that  using  Q  =  P  +  2  worked  extremely  well  for  all  of  the  test  cases  using  any 
grid.  Also,  for  P  >  4  the  accuracy  differences  are  insignificant  between  inverting  the  mass 
matrix  in  its  full  form  or  in  its  diagonal  form.  For  this  reason,  the  mass  matrix  on  the  left-hand 
side  is  stored  in  its  diagonal  form.  This  is  what  is  normally  done  in  spectral  element  methods; 
the  justification  is  that  Q  =  P  is  normally  used  which  then  yields  a  diagonal  matrix  anyway. 
We  are  currently  using  Q  =  P  +  2  because  it  appears  to  be  a  bit  more  stable  especially  for  the 
longer  runs.  Using  Q  =  P  requires  the  use  of  filters  (see  Reference  [6]);  in  our  runs  with 
Q  =  P  +2  no  filtering  of  any  kind  was  required. 


3.2.  Time  discretization 

Starting  from  the  spatial  discretization  given  in  Equation  (4),  the  resulting  system  of  ordinary 
differential  equations  (ODEs)  can  now  be  written  as 


dip 

dt 


H(<p) 


where 


H(<p)  =  Mq  *[(A ijk-ak)(pj  +  St(ip )] 

Recall  that  the  matrix  M  is  a  diagonal  matrix  and  so  it  requires  very  little  storage  and  does  not 
require  the  inversion  of  a  matrix.  After  integrating  in  time  by  a  general  family  of  Runge-Kutta 
schemes  yields 
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<pk+1  =  <p"  +  At/3H(<pf-1 
where 

/?  =  7  ,  ,  ,  k=\,...,K  and  H((p)°  =  H(<p)n 

K  —  k  +  1 

This  one-step  multi-stage  method  is  used  for  the  first  couple  of  iterations,  at  which  point  the 
more  efficient  (but  slightly  less  accurate)  third-order  Adams-Bashforth  scheme  is  used.  The 
time  discretization  then  is 

<pn  +  1  =  <p"  +  ^[23H((p)n-  16H(<p)"-1  +  5H(<p)"-2] 

which  clearly  requires  knowing  the  conservation  variables  tp  at  the  three  time  levels  n,  n  —  1 
and  n  —  2.  For  this  reason,  we  use  the  Runge-Kutta  scheme  for  the  first  two  time  steps. 

The  third-order  Adams-Bashforth  scheme  has  some  advantages  and  disadvantages  over 
other  explicit  schemes,  such  as  the  third-  or  fourth-order  Runge-Kutta  schemes.  The  main 
drawback  of  the  Adams-Bashforth  scheme  is  that  it  requires  a  very  small  time  step  for 
stability  purposes  as  compared  with  the  Runge-Kutta  scheme;  however,  the  main  advantage 
is  that  this  scheme  is  much  faster  than  the  Runge-Kutta  scheme.  By  this  reason  alone,  the 
Adams-Bashforth  scheme  is  usually  selected  as  the  time  integration  scheme  of  choice  for 
spectral  element  models  [2,5,6]. 

Iskandarani  [2]  and  Ma  [5]  and  have  shown  that  the  stability  of  the  scheme  requires  the  time 
step  limit  to  be 

(14) 

where  p  is  the  order  of  the  polynomial,  Ne  is  the  number  of  elements,  and  C  is  a  proportion¬ 
ality  constant.  Therefore,  for  a  fixed  number  of  elements,  increasing  the  order  of  the 
polynomial  by  a  factor  of  2  requires  the  time  step  to  be  decreased  by  a  factor  of  4.  This 
relation  held  for  all  our  numerical  experiments  and  in  fact  was  used  to  determine  the  time  step 
size  for  a  given  p  and  Nc. 

3.2.1.  Lagrange  multiplier  constraint.  Once  the  solution  has  been  obtained  at  the  new  time  step, 
we  need  to  constrain  the  fluid  particles  to  remain  on  the  surface  of  the  sphere.  We  then  have 
for  the  momentum  equations 

(^u)2+1  =(<pu)u+1  +  pr 

where  the  subscripts  ‘c’  and  V  denote  the  constrained  and  unconstrained  values.  Taking  the 
scalar  product  of  r,  where  r  is  the  co-ordinate  vector  of  each  grid  point,  gives 
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r-((pu)nc  +  l  =  r-(<pu)"+1  +  fia2  (15) 

For  a  fluid  particle  to  remain  on  the  surface  of  the  sphere,  the  velocity  field  must  be 
orthogonal  to  the  co-ordinate  (radial)  vector  of  its  grid  point.  In  other  words 

u-r  =  0 

This  means  that  the  left-hand  side  term  in  Equation  (15)  vanishes,  and  after  rearranging,  yields 
r‘(<pu)”+ 1 


4.  ICOSAHEDRAL  GRIDS 

One  of  the  difficulties  of  using  the  typical  icosahedral  grid  used  in  References  [7,13]  is  that  it 
is  restricted  to  only  a  few  grids  because  the  number  of  grid  points  increases  by  a  factor  of  4 
with  each  refinement  of  the  grid.  In  other  words,  from  q  =  4  to  5  (where  q  is  the  refinement 
level)  we  go  from  2562  grid  points  to  10242  grid  points.  This  is  a  very  limiting  constraint  when 
trying  to  do,  say,  a  convergence  study  on  a  computer  of  modest  size.  One  way  around  this 
dilemma  is  to  construct  the  general  icosahedral  grid  we  developed  in  Reference  [12]. 

We  begin  with  the  initial  icosahedron  having  12  points  and  20  equilateral  triangular 
elements.  Then  we  subdivide  each  triangle  of  the  initial  icosahedron  by  an  nth  order  Lagrange 
polynomial.  Before  doing  so,  however,  it  is  best  to  map  this  triangle  onto  a  gnomonic  space. 
The  most  unbiased  mapping  is  obtained  by  mapping  about  the  centroid  of  the  icosahedral 
triangle.  (Note  that  because  all  of  the  initial  triangles  are  equilateral  then  the  centroid  is 
equivalent  to  the  circumcenter  of  the  triangle.)  Let  Uc,  6C)  be  the  centroid  of  the  triangle  we 
wish  to  map.  Then,  the  gnomonic  mapping  is  given  by 

a  cos  6  sin(2  —  Ac)  «[cos  0C  sin  6  —  sin  8C  cos  6  cos(2  —  2C)] 

sin  6C  sin  0  +  cos  8C  cos  0  cos(/l  —  Ac )  '  sin  6C  sin  6  +  cos  8C  cos  0  cos(/l  —  Ac ) 

(16) 

which  is  rather  complicated.  However,  by  first  applying  a  rotation  mapping  whereby  in  the 
new  co-ordinate  system  the  co-ordinates  (/,  0)  are  located  at  (0,0),  then  Equation  (16) 
becomes 

x  =  atmA',  y  =  a  tan  6'  sec  A'  (17) 

where  the  rotation  mapping  is  given  by 
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cos  0  sm(x  —  ic) 

/  =  arctan  — — - — ; — - - - — - - - — 

_sm  0C  sin  u  +  cos  0C  cos  0  cos(/  —  Ac)_ 

9'  =  arcsin[cos  9C  sin  0  —  sin  9C  cos  0  cos(x  —  2C)]  (18) 

This  approach  results  in  the  construction  of  a  general  icosahedral  grid  defined  by 
Nj  =  10(77  -  l)2  +  20(71  -  1)  +  12 
N]  =  2(Nl-2) 

Nj  =  3(NTp-2) 

where  Np,  Nj,  and  Nj  denote  the  number  of  points,  elements,  and  sides  comprising  the 
triangular  grid,  and  n  is  the  order  of  the  Lagrange  polynomial  used  to  subdivide  the  20  initial 
triangles  of  the  icosahedron.  This  grid  [12]  was  shown  to  be  much  more  general  and  allows 
much  more  flexibility  in  constructing  grids  than  the  typical  icosahedral  grid  used  in  References 
[7,13],  Note  that  if  we  write  n  =  2q  we  recover  the  typical  icosahedral  grid. 

Once  we  have  constructed  the  triangular  icosahedral  grid,  we  then  subdivide  each  triangular 
element  into  three  smaller  quadrilateral  elements.  We  need  quadrilaterals  because  the  high- 
order  nodal  basis  functions  are  only  defined  on  quadrilaterals  and  not  on  triangles  (\j/  in 
Equation  (9)).  Upon  dividing  the  triangles  into  quadrilaterals  we  can  then  construct  the 
spectral  collocation  points  inside  each  element  resulting  in  a  quadrilateral  grid  with  the 
following  properties: 

Af  =  6(A'2-2)(p)2  +  2 

Nf  =  6(7Vp  —  2) 

Af  =  12(Ap'-2) 

where  N Nf,  and  Nf  denote  the  number  of  points,  elements,  and  sides  comprising  the 
quadrilateral  spectral  element  grid,  and  p  is  the  order  of  the  spectral  element  basis  functions. 
Table  I  gives  the  grid  dimensions  for  various  values  of  n  and  p.  Figure  1  shows  some  possible 
grid  configurations  for  n  =  1  and  p  =  1, .  . .  ,  32. 

This,  of  course,  is  not  the  only  possible  grid  that  can  be  used  for  spectral  elements  on  the 
sphere.  Alternatively,  we  could  use  the  cubic  gnomonic  grids  employed  in  References  [6,21].  In 
fact,  this  grid  appears  to  be  a  much  more  natural  choice  because  the  grid  is  constructed  by 
mapping  a  hexahedron  (cube)  onto  the  sphere.  Each  of  the  six  faces  can  then  be  subdivided 
into  as  many  quadrilaterals  as  required.  The  main  drawback  of  this  grid  is  that  the 
quadrilateral  elements  furthest  from  the  centroid  of  each  face  are  much  smaller  than  those  near 
the  centroid  due  to  the  distortion  of  the  gnomonic  projection  given  by  Equation  (16).  Even 
though  we  are  also  using  this  projection  for  our  icosahedral  grid  we  now  have  20  smaller 
elements  to  map  rather  than  six  larger  ones  for  the  cubic  gnomonic  grid.  To  illustrate  the 
point,  comparing  both  grids  for  a  comparable  number  of  elements  we  find  that  the  ratio  of 
maximum  to  minimum  element  sizes  ranges  from  1  to  1.96  for  the  icosahedral  grid  and  from 
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Table  I.  The  number  of  grid  points  (N®),  elements  (/VJ?),  and  sides  (N$)  for 
the  spectral  element  icosahedral  grid  using  quadrilaterals. 


n 

P 

n  *  p 

nq 

v  p 

N? 

N ? 

4 

1 

4 

962 

960 

1920 

8 

1 

8 

3842 

3840 

7680 

16 

1 

16 

15  362 

15  360 

30  720 

32 

1 

32 

61  442 

61  440 

122  880 

64 

1 

64 

245  762 

245  760 

491  520 

1 

4 

4 

962 

60 

120 

2 

4 

8 

3842 

240 

480 

4 

4 

16 

15  362 

960 

1920 

8 

4 

32 

61  442 

3840 

7680 

16 

4 

64 

245  762 

15  360 

30  720 

1 

8 

8 

3842 

60 

120 

2 

8 

16 

15  362 

240 

480 

4 

8 

32 

61  442 

960 

1920 

8 

8 

64 

245  762 

3840 

7680 

1 

16 

16 

15  362 

60 

120 

2 

16 

32 

61  442 

240 

480 

4 

16 

64 

245  762 

960 

1920 

1 

4 

4 

962 

60 

120 

1 

8 

8 

3842 

60 

120 

1 

16 

16 

15  362 

60 

120 

1 

32 

32 

61  442 

60 

120 

1 

64 

64 

245  762 

60 

120 

The  order  of  the  generalized  icosahedral  grid  is  given  by  n  and  the  order  of  the  spectral 
element  basis  function  is  given  by  p. 


2.20  to  5.14  for  the  cubic  grid  (corresponding  to  n  =  1, .  .  . ,  64).  Clearly,  the  icosahedral  grid 
yields  a  more  uniform  representation  on  the  sphere  which  is  a  highly  desirable  property  for 
global  atmospheric  models;  the  construction  of  a  global  atmospheric  model  is  the  long  term 
goal  of  our  research. 


5.  RESULTS 

For  the  numerical  experiments,  the  following  three  terms  are  used  for  judging  the  performance 
of  the  scheme: 

the  L2  error  norm 


Published  in  2001  by  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Fluids  2001;  35:  869-901 


» 


»I»] 

*1*1 


884 


F.  X.  GIRALDO 


where  <p  represents  the  conservation  variables,  and  the  following  two  additional  conservation 
measures: 


M  = 


- 

cp  dQ 

7s  > 


exact 

Jn 


E  = 


<p(u2  +  v2  +  w2)  +  cp2  dQ 

Jq 

r* 

^ exactly  exact  1  ^  exact  4”  W  exact)  b  exact 

Jn 


The  L2  error  norm  compares  the  root-mean-square  (r.m.s.)  per  cent  error  of  the  numerical  and 
exact  solutions,  M  measures  the  conservation  property  of  the  mass,  and  E  measures  the 
conservation  of  the  total  available  energy.  The  ideal  scheme  should  yield  an  L2  error  norm  of 
zero,  and  mass  and  energy  measures  of  one. 

Six  test  cases  are  used  in  order  to  test  the  algorithm.  Test  cases  1,  2,  3,  5,  and  6  correspond 
to  the  test  cases  given  in  Reference  [22].  Test  case  4  is  given  in  Reference  [23],  Test  case  1 
involves  the  mass  equation  only,  whereas  the  remainder  of  the  test  cases  concern  the  full 
shallow  water  equations.  In  addition,  cases  1-3  have  analytic  solutions  and  are  used  to 
determine  the  accuracy  of  the  spectral  element  method  quantitatively.  Test  cases  4-6,  on  the 
other  hand,  do  not  have  analytic  solutions  and  are  thus  only  used  to  determine  the  accuracy 
of  the  scheme  qualitatively. 

For  all  of  the  test  cases  studied  in  this  paper,  the  time  steps  used  were  A?  =  864,  216,  54, 
and  14  s  for  the  grids  (n  =  1)  p  =  4,  8,  16,  and  32  respectively.  Note  that  the  time  step  size 
decreases  by  a  factor  of  4  for  p  increasing  from  2'"  to  2'” + 1  for  m  =  2, .  .  . ,  4  as  given  by  the 
Adams- Bashforth  stability  limit  (Equation  (14)). 


5.1.  Case  1:  steady  state  advection 

This  test  case  concerns  the  solid  body  rotation  of  a  cosine  wave.  It  only  tests  the  mass  equation 
as  the  velocity  field  is  assumed  to  remain  unchanged  throughout  the  computation.  Results  are 
obtained  after  one  full  revolution,  which  corresponds  to  an  integration  of  12  days. 

By  using  the  mapping  from  spherical  to  Cartesian  space 


x  =  a  cos  9  cos  A 


y  =  a  cos  0  sin  A 
z  =  a  sin  0 
where 


A  =  arctam 


(*) 


we  can  write  the  initial  conditions  in  terms  of  Cartesian  co-ordinates.  The  spherical  co¬ 
ordinates  (A,  9)  correspond  to  the  longitude  and  latitude.  This  results  in  the  following  velocity 
field  in  Cartesian  space 
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u  =  —  us  sin  X  —  vs  sin  6  cos  X 
v  =  +  us  cos  X  —  vs  sin  6  sin  X 
w  =  +  vs  cos  6 

where  us  and  vs  are  the  zonal  and  meridional  velocity  components  in  spherical  co-ordinates. 

A  convergence  study  is  shown  in  Figure  2  for  the  schemes  p  =  1,  p  =  4,  p  =  8,  p  =  16,  and 
n  =  1 .  The  L2  error  of  the  mass  is  plotted  against  the  product  n  *p.  The  scheme  p  =  1  represents 
keeping  the  order  of  the  polynomial  fixed  at  first-order  while  increasing  the  number  of 
elements  (increasing  n  from  4  to  32).  Likewise,  p  =  4  represents  keeping  the  order  of  the 
polynomial  fixed  at  fourth  order  while  increasing  the  number  of  elements  (increasing  n  from 
1  to  8).  Similarly,  p  =  8  represents  keeping  the  order  of  the  polynomial  fixed  at  eighth  order 
while  increasing  the  number  of  elements  (increasing  n  from  1  to  4),  and  p=  16  represents 
keeping  the  order  of  the  polynomial  fixed  at  sixteenth  order  while  increasing  the  number  of 
elements  (increasing  n  from  1  to  2).  Finally,  n  =  1  represents  keeping  the  number  of  elements 
fixed  (n  =  1)  but  increasing  the  order  of  the  polynomial  (increasing  p  from  4  to  32).  Table  I 
contains  the  number  of  grid  points  and  elements  for  all  of  the  grid  combinations  of  p  and  n 
used  throughout  the  paper. 


Figure  2.  Case  1 .  Convergence  study  of  the  L2  mass  error  for  various  values  of  n  *p. 
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This  study  illustrates  that  a  better  convergence  rate  is  achieved  by  keeping  the  number  of 
elements  fixed  and  increasing  the  order  of  the  polynomial  (n  =  1).  The  p  =  1  scheme  converges 
at  a  rate  of  2.15,  meaning  that  by  doubling  the  variable  n  *p  we  decrease  the  error  by  a  factor 
of  2.15,  while  the  p  =  4,  p  =  8,  p  =  16,  and  n  =  1  schemes  converge  at  rates  of  5.85,  6.32,  6.60, 
and  10  respectively.  The  n  =  1  scheme  convergences  at  a  spectral  (exponential)  rate;  this  is  the 
p-version  of  the  spectral  element  method.  The  p  =  1  scheme  gives  the  slowest  convergence  rate 
and  it  represents  using  linear  elements  but  increasing  the  number  of  elements  (increasing  n 
from  4  to  32);  this  is  the  /(-version  of  the  spectral  element  method  and,  in  fact,  it  is  identical 
to  a  bilinear  finite  element  method.  For  smooth  flows  (which  this  test  case  represents),  the 
/(-version  is  superior  to  the  /? -version  but  for  general  problems  a  combination  of  both  (the  h  p 
method)  is  optimal.  For  this  reason  we  have  also  included  the  p  =  4,  p  =  8,  and  p  =  16  schemes. 


5.2.  Case  2:  global  steady  state  non-linear  zonal  geostrophic  flow 

This  case  is  a  steady  state  solution  to  the  non-linear  shallow  water  equations.  The  equations 
are  geostrophically  balanced  and  remain  so  for  the  duration  of  the  integration.  The  velocity 
field  thus  remains  constant  throughout  the  computation.  The  geopotential  height  (q>)  under¬ 
goes  a  solid  body  rotation,  but  since  the  initial  height  field  is  given  as  a  constant  band,  the 
solution  looks  the  same  throughout  the  time  integration.  The  velocity  field  is  the  same  as  that 
used  in  case  1.  All  the  results  reported  are  for  a  5-day  integration  as  suggested  in  Reference 
[22]. 

A  convergence  study  is  shown  in  Figure  3  for  the  p=  1,  p  =  4,  and  n  =  1  schemes.  The 
convergence  rates  for  the  three  schemes  p=  1,  p  =  4,  and  n=  1,  are  3.73,  6.50,  and  8.74 
respectively.  Once  again,  the  p -version  of  the  spectral  element  method  (n  =  1)  yields  superior 
results.  The  n  =  1  scheme  falls  slightly  short  of  exponential  convergence  but  it  gives  impressive 
results  nonetheless.  However,  it  would  appear  that  the  scheme  reaches  machine  precision  and 
hence  plateaus  for  increasing  n*p  values.  For  this  reason,  the  p  =  8  and  p  =  16  schemes  were 
omitted  from  this  plot  because  the  results  were  virtually  indistinguishable  from  the  n  =  1 
results.  The  convergence  rate  of  the  n  =  1  scheme  is  higher  than  exponential  for  p  <  8  but  slows 
down  after  8.  In  fact,  from  p  =  4  to  p  =  8  the  scheme  is  superconvergent  with  the  error 
decreasing  by  three  orders  of  magnitude! 

5.3.  Case  3:  steady  state  non-linear  zonal  geostrophic  flow  with  compact  support 

This  case  is  similar  to  case  2  except  that  the  velocity  is  zero  everywhere  except  in  a  very  small 
isolated  region.  This  isolated  region,  or  jet,  encapsulates  the  flow  and  limits  the  geopotential 
height  field  to  remain  within  a  very  confined  circular  region.  As  in  case  2,  the  results 
correspond  to  a  5-day  integration. 

A  convergence  study  is  shown  in  Figure  4  for  the  p  =  1,  p  =  4,  p  =  8,  p  =  16,  and  n  =  1 
schemes.  The  convergence  rates  for  the  five  schemes  are  2  for  the  p  =  1,  p  =  4,  p  =  8,  and 
p  =  16  schemes,  while  it  is  3.59  for  the  n  =  1  scheme.  Although  the  convergence  rates  for  this 
test  case  are  low  relative  to  the  previous  two  test  cases  it  does  show  that  even  for  localized 
phenomena  (such  as  this  test  case  represents)  the  n  =  \  scheme  yields  the  best  results, 
converging  at  the  fastest  rate. 
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Figure  3.  Case  2.  Convergence  study  of  the  L2  mass  error  for  various  values  of  n  *p. 

Having  performed  satisfactorily  in  these  three  analytic  test  cases,  we  now  move  to  three 
more  test  cases  without  analytic  solutions.  These  test  cases  shall  be  used  to  gauge  the 
qualitative  performance  of  the  spectral  element  method. 

5.4.  Case  4:  dancing  high-low  waves 

This  test  case  comes  from  Reference  [23]  and  is  not  an  analytic  solution  to  the  shallow  water 
equations.  The  initial  geopotential  height  is  comprised  of  two  large  waves  with  the  left  wave 
being  the  low  wave  and  the  right  wave  being  the  high  wave,  when  viewed  from  the  north  pole 
(see  Figures  5-7).  The  waves  rotate  clockwise  in  a  dance -fashion  so  that  after  5  days  of 
integration,  the  high  wave  is  now  on  the  left  and  the  low  wave  is  on  the  right. 

Figures  5-7  show  the  mass  and  velocity  field  after  5  days  for  the  three  grid  resolutions  of 
p  =  4,  p  =  8,  and  p=  16  (with  n=  1  for  all  three  grids).  The  viewpoint  in  these  figures  is 
(A,  9)  =  (0,  90),  which  represents  a  view  from  above  the  north  pole.  We  can  see  from  these 
three  figures  that  the  contours  have  a  similar  structure  for  the  three  grids  but  the  flow  features 
become  better  resolved  as  p  is  increased.  For  instance,  comparing  the  mass  contour  plots  for 
p  =  4  and  p  =  8  (Figures  5(a)  and  6(a))  we  can  see  that  the  contour  curves  have  become  sharper 
everywhere.  The  resolution  further  increases  as  we  go  from  =  8  to  p  =  16  (Figures  6(a)  and 
7(a)).  The  increase  in  resolution  is  even  more  apparent  with  the  zonal  and  meridional  velocity 
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Figure  4.  Case  3.  Convergence  study  of  the  L2  mass  error  for  various  values  of  n  *p. 


Figure  5.  Case  4.  The  (a)  mass,  (b)  zonal  velocity,  and  (c)  meridional  velocity  on  grid  n  =  1,  p  =  4  after 

5  days  viewed  from  (i,  8)  =  (0,  90). 

components.  Figure  7(b)  and  (c)  show  a  much  sharper  wave  structure  than  those  in  Figure  5(b) 
and  (c)  and  Figure  6(b)  and  (c).  The  difference  in  resolution  is  more  apparent  near  the  center 
of  these  plots  (the  north  pole)  because  this  is  the  center  of  rotation  of  the  two  waves;  large 
gradients  exist  at  the  center  where  the  low  and  high  waves  meet.  Clearly,  this  test  case  shows 


Published  in  2001  by  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Fluids  2001;  35:  869-901 


A  SPECTRAL  ELEMENT  SHALLOW  WATER  MODEL  ON  SPHERICAL  GEODESIC  GRIDS 


889 


a) 


b) 


c) 


Figure  6.  Case  4.  The  (a)  mass,  (b)  zonal  velocity,  and  (c)  meridional  velocity  on  grid  n  =  1,  p  =  8  after 

5  days  viewed  from  (2,  8)  =  (0,  90). 


Figure  7.  Case  4.  The  (a)  mass,  (b)  zonal  velocity,  and  (c)  meridional  velocity  on  grid  n  =  1,  p  =  16  after 

5  days  viewed  from  (2,  9)  =  (0,  90). 

that  the  spectral  element  method  is  converging  towards  a  wave  structure  best  represented  by 
Figure  7. 

We  have  run  this  test  case  for  a  1 5-day  integration  and  the  waves  remain  intact  even  at  this 
long  integration  period.  However,  we  only  show  results  after  5  days  because  the  reference 
solution  in  Reference  [22]  only  shows  results  for  a  5-day  integration.  Furthermore,  the  results 
illustrated  in  Figures  5(a),  6(a),  and  7(a)  are  quite  similar  to  Figures  5  and  6  in  Reference  [23], 
The  results  presented  in  Reference  [23]  were  obtained  by  a  finite  difference  semi-Lagrangian 
method  using  7200  points.  Because  our  spectral  element  method  is  higher  order  and  we  used 
15362  points  (for  p  =  16),  differences  between  the  two  solutions  are  expected. 

Because  the  above  contour  plots  only  show  a  hemisphere,  we  have  plotted  the  results  of  the 
p  =  16  grid  on  a  longitude-latitude  projection  in  Figure  8,  thereby  showing  the  resulting  wave 
structure  throughout  the  entire  globe.  The  longitude-latitude  contour  plot  is  shown  for  a 
5-day  integration.  Note  that  the  longitudinal  view  in  this  plot  is  2  =  180,  whereas  in  Figures 
5-7  it  is  2  =  0.  For  this  reason,  the  high  and  low  waves  are  reversed  in  Figure  8.  In  the 
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Figure  8.  Case  4.  Contour  plot  of  the  mass  for  the  grid  n  =  1,  p  =  16  after  5  days  on  a  longitude-latitude 

(/,  6)  projection. 

northern  hemisphere  of  this  figure,  the  high  wave  is  clearly  seen  on  the  right  and  the  low  wave 
on  the  left  and  in  the  southern  hemisphere  their  positions  are  reversed.  Note  that  the  high 
waves  in  the  northern  and  southern  hemisphere  have  the  same  values.  The  same  holds  true  for 
the  low  waves.  This  implies  that  the  numerical  method  and  grid  are  maintaining  the  initial 
symmetries  of  the  problem  throughout  the  5-day  integration  without  introducing  any  numeri¬ 
cal  biasing  to  either  hemisphere. 

5.5.  Case  5:  zonal  flow  over  an  isolated  mountain 

This  case  is  similar  to  case  2  except  that  now  there  is  a  mountain  on  the  sphere.  This  is  the 
only  problem  in  the  test  cases  studied  here  which  includes  topography.  The  mountain  is  conical 
in  shape  and  forces  the  zonal  flow  to  deflect  off  the  mountain.  Due  to  the  zonal  flow 
impinging  on  the  mountain,  various  wave  structures  form  and  propagate  around  the  sphere. 
The  results  shown  for  this  test  case  are  for  1-,  5-,  and  10-day  integrations. 

Figures  9-11  show  the  convergence  of  the  spectral  element  method  for  the  three  different 
resolutions  p  =  4,  p  =  8,  and  p  =  16  (with  n  =  1).  The  viewpoint  used  for  these  plots  is 
(A,  0)  =  (270,  0).  This  viewpoint  was  chosen  because  the  mountain  is  located  at  (A,  6)  = 
(270,  30)  (see  Reference  [22]  for  details  of  this  test  case).  At  1  day,  the  wave  structure  generated 
by  the  mountain  is  immediately  noticeable  in  Figures  9(a),  10(a),  and  11(a).  However,  for 
p=  16  (Figure  11(a))  the  wave  structure  is  much  sharper  than  for  the  other  two  resolutions. 
For  p  =  4  we  get  some  semblance  of  the  wave  structure  (Figure  9)  but  note  that  the  mass 
contour  plots  become  more  jagged  at  5  and  10  days.  Figure  10  shows  that  for  p  =  8  the 
contours  are  smoother,  but  they  nonetheless  lose  some  smoothness  near  the  north  pole  and  at 
the  equatorial  troughs,  especially  after  10  days.  In  contrast,  the  sharpness  in  resolution 
achieved  by  the  highest-order  scheme,;?  =  16  (Figure  11),  is  evident  at  both  5  and  10  days.  For 
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Figure  9.  Case  5.  The  mass  at  (a)  1,  (b)  5,  and  (c)  10  days  on  grid  n=  1,  p  =  4  viewed  from 

(A,  0)  =  (270,  0). 


Figure  10.  Case  5.  The  mass  at  (a)  1,  (b)  5,  and  (c)  10  days  on  grid  n=  1,  p  =  8  viewed  from 

(A,  6)  =  (210,  0). 


Figure  11.  Case  5.  The  mass  at  (a)  1,  (b)  5,  and  (c)  10  days  on  grid  n=  1,  p=  16  viewed  from 

(A,  6)  =  (210,  0). 
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example,  in  the  vicinity  of  the  mountain,  we  observe  the  contours  to  be  quite  sharp.  In  addition, 
the  equatorial  troughs  are  well  defined  and  the  north  polar  circulation  triggered  by  the  mountain 
is  quite  smooth.  It  should  be  mentioned  that  no  type  of  smoothing  has  been  used  for  any  of 
the  results  illustrated — either  for  plotting  or  for  damping  the  high-order  wave  numbers  of  the 
numerical  scheme.  The  smoothness  observed  in  the  contour  plots  is  a  result  of  the  increased  grid 
resolution;  in  these  cases,  obtained  by  increasing  p  only  (the  polynomial  order  of  the  spectral 
element  basis  functions). 

Figure  12  shows  the  contour  plot  of  the  mass  on  a  longitude-latitude  projection  for  p  =  16 
for  a  10-day  integration.  These  results  show  the  same  wave  structure  given  in  References  [7,24]; 
it  should  be  mentioned  that  the  current  results  have  been  shifted  180°  in  order  to  compare  against 
their  results  because  in  the  current  case  the  mountain  is  located  at  X  =  270  (as  suggested  in 
Reference  [22]),  whereas  in  theirs  it  is  stationed  at  X  =  90.  In  fact,  the  contours  in  Figure  12  match 
those  in  figure  5.  lc  of  Reference  [24]  exactly.  The  results  presented  in  Reference  [24]  were  obtained 
using  a  spectral  method  with  a  resolution  of  T213,  consisting  of  204800  points.  Our  results  were 
obtained  with  15362  points  (p  =  16  with  n  =  1). 

5.6.  Case  6:  Rossby-Haurwitz  wave 

Although  Rossby-Haurwitz  waves  are  not  analytic  solutions  to  the  shallow  water  equations, 
they  have  become  a  de  facto  test  case.  In  a  non-divergent  barotropic  model,  the  initial 
geopotential  height  field  undergoes  a  solid  body  rotation  in  a  counterclockwise  direction  (when 
viewed  from  the  north  pole).  The  angular  velocity  is  given  by 

_  R( 3  +  R)co  -  2Q 
(1  -  R )(2  •  A') 

where  R  =  4  is  the  wavenumber.  Results  are  reported  for  1-  and  7-day  integrations. 


Figure  12.  Case  5.  Contour  plot  of  the  mass  for  the  grid  n  =  1,  p  =  16  after  10  days  on  a  longitude- 

latitude  (X,  6)  projection. 
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Figures  13-15  show  the  mass  contours  for  the  three  different  grids  from  the  viewpoint 
(A,  9)  =  (0,  90).  This  particular  viewpoint  was  selected  because  from  this  vantage  point,  the 
waves  undergo  a  counterclockwise  solid  body  rotation,  which  then  simplifies  comparisons 
between  different  resolution  grids.  After  1  day,  the  wave  is  completely  intact  for  all  three 
resolutions,  but  as  the  resolution  is  increased  a  sharper  wave  image  is  obtained.  At  7  days,  the 
p  =  4  grid  begins  to  dissolve  the  wave  structure,  whereas  it  remains  completely  intact  for  the 
p  =  8  and  p  =  16  grids.  In  addition,  the  p  =  16  grid  gives  a  sharper  image  than  the  p  =  8  grid. 

In  Figure  16  we  have  plotted  the  mass  contours  on  a  longitude -latitude  projection  for  the 
grid  p  =  16  after  7  days  in  order  to  show  the  full  wave  structure  throughout  the  sphere.  The 
results  in  Figure  16  match  those  in  figure  5.5c  of  Reference  [24]  exactly.  Once  again,  the  results 
given  in  Reference  [24]  were  obtained  using  a  spectral  method  with  a  resolution  of  T213,  which 
is  comprised  of  204800  points.  In  comparison,  our  p  =  16  grid  (with  n  =  1)  contains  15362 
points.  Note  that  the  solution  after  7  days  is  completely  symmetrical,  thereby  confirming  that 
the  scheme  has  conserved  the  symmetry  of  the  initial  conditions. 

5.7.  Summary  of  results 

The  spectral  element  method  on  icosahedral  grids  performed  extremely  well  in  all  six  test  cases. 
By  discretizing  the  conservation  form  of  the  equations  with  the  spectral  element  method  we 
were  able  to  conserve  mass  (M)  exactly  for  all  six  test  cases.  The  total  energy  (E)  was 
conserved  exactly  for  cases  1-3.  For  cases  4-6  the  spectral  element  method  lost  some  of  its 
initial  energy,  with  the  worst  case  representing  only  a  very  small  loss  (0.1  per  cent  for  a  15-day 
integration). 


x 


Figure  16.  Case  6.  Contour  plot  of  the  mass  for  the  grid  n=  1,  p=  16  after  7  days  on  a  longitude- 

latitude  (2,  9)  projection. 
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The  results  for  all  six  test  cases  show  that  the  spectral  element  discretization  of  the  Cartesian 
conservation  form  of  the  shallow  water  equations  on  icosahedral  grids  is  extremely  accurate. 
However,  these  results  only  prove  without  a  doubt  that  the  method  is  stable  for  at  least  15 
days.  Case  1  is  clearly  stable  for  any  length  of  time  because  it  only  involves  the  mass  equation. 
Cases  2-6  were  run  up  to  15  days  without  any  instabilities  arising. 

Currently,  the  numerical  algorithm  does  not  use  any  type  of  filtering.  Spectral  element 
models  typically  use  some  kind  of  filtering  because  this  method,  like  the  spectral  method,  is  not 
immune  to  aliasing  errors  that  arise  from  the  nonlinear  terms.  It  is  possible  that  the 
higher-order  schemes  (p  >  8)  of  our  spectral  element  method  would  require  some  form  of 
filtering  in  order  to  suppress  the  growth  of  spurious  modes  for  very  long  time  integrations. 
Further  testing  of  the  model  should  reveal  whether  or  not  filtering  is  required.  If  so,  then  we 
expect  to  explore  ad  hoc  spectral  element  filters  like  those  presented  in  References  [6,25]. 

5.8.  Computational  cost 

Figure  17  summarizes  the  computational  cost  of  using  the  spectral  element  method  on 
icosahedral  grids  for  the  p  =  1,  p  =  4,  p  =  8,  p  =  16,  and  n  =  1  schemes.  These  results  were 
obtained  using  case  3  but  are  generally  representative  of  our  spectral  element  method.  These 


Figure  17.  The  total  computational  time  incurred  by  the  spectral  element  method  as  a  function  of  grid 
parameter  n*p  for  case  3  for  a  5-day  integration. 
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results  were  obtained  on  a  DEC  Alpha  8400  (EV5  CPU  chip)  with  a  clock  speed  of  300  MHz 
running  six  threads  (multi-threading)  using  Open  MP. 

It  is  obvious  from  looking  at  this  figure  that  the  A -version  of  our  method  (p  =  1)  costs  the 
least  while  the  p-version  (n  =  1)  costs  the  most.  The  reason  for  this  is  due  to  the  spectral 
element  method  requiring  0(n2p3)  operations  per  time  step.  Since  in  the  A -met  hod  we  increase 
n  while  keeping  p  fixed,  the  number  of  operations  increases  by  powers  of  2,  while  in  the 
p -method  we  increase  p  and  hold  n  constant  thereby  increasing  the  number  of  operations  by 
powers  of  3. 

This  figure  shows  that  the  cost  of  running  the  p  =  4  and  p  =  1  methods  are  almost  identical. 
In  fact,  increasing  p  from  1  to  8  only  increases  the  cost  by  a  factor  of  1.5.  Increasing;;  to  16 
increases  the  cost  by  2.5.  Finally,  using  the  n  =  1  method  clearly  costs  much  more,  in  fact 
almost  an  order  of  magnitude  more.  These  numbers  are  for  n*p  =  32  in  Figure  17. 

So  the  issue  is  to  determine  what  is  the  best  choice  of  n  and  p  which  yields  the  best  results 
while  doing  so  in  the  most  efficient  manner  possible.  Revisiting  Figure  4  we  see  (again  for 
n  *p  =  32)  that  going  from  p  =  1  to  p  =  8  decreases  the  error  by  a  factor  of  3  while  only 
increasing  the  cost  by  a  factor  of  1.5.  Increasing;;  to  16  reduces  the  error  by  a  factor  of  5  but 
only  increasing  the  cost  by  a  factor  of  2.5.  Using  the  n  =  1  scheme  (i.e.,  ;;  =  32)  decreases  the 
error  by  an  order  of  magnitude.  Of  course  in  this  case  the  cost  has  also  increased  by  an  order 
of  magnitude.  Therefore,  the  optimal  range  of ;;  must  be  somewhere  between  8  and  32  but  we 
can  take  the  more  conservative  range  of  8-16. 

Looking  at  it  another  way,  we  can  argue  quite  convincingly  that  the  high-order  method 
(p  >  8)  is  much  more  efficient  than  the  low-order  method  (;;  =  1).  In  Figure  4  the  ;;  value  for 
the  n  =  1  curve  required  to  yield  the  same  error  as  that  for  the  ;;  =  1  curve  for  n*p  =  32  is 
approximately  p  =  10  (for  an  error  of  10  ~3).  From  Figure  17  the  n=  1  scheme  with  ;;  =  10 
costs  an  order  of  magnitude  less  than  the  ;;  =  1  scheme  with  n  =  32  (2  x  102  versus  4  x  1  O’  s). 
Therefore,  the  high-order  methods  give  not  only  increased  accuracy  but  increased  efficiency  as 
well. 

The  algorithm  currently  does  not  scale  linearly  with  increasing  ;;.  Certainly  for  ;;  =  32  the 
method  returns  as  much  accuracy  as  it  costs.  Going  beyond  this  value  will  incur  even  higher 
costs.  However,  larger  values  of  ;;  (  >  32)  may  become  practical  if  the  algorithm  can  be 
streamlined  and  massively  parallel  architectures  exploited.  We  are  currently  running  an  Open 
MP  (multi-threading)  version  of  the  model  but  we  hope  to  implement  a  message  passing 
interface  (MPI)  version  in  the  near  future.  In  fact,  due  to  the  explicit  time  differencing 
currently  used  it  should  be  trivial  to  develop  an  MPI  version. 


6.  CONCLUSIONS 

The  results  obtained  for  all  six  problems  show  that  the  spectral  element  method  on  icosahedral 
grids  offers  a  viable  method  for  solving  the  spherical  shallow  water  equations  in  Cartesian 
conservation  form.  Rather  than  solving  the  equations  in  spherical  co-ordinates,  the  equations 
are  transformed  to  Cartesian  co-ordinates,  which  have  two  advantages  over  spherical  co-ordi¬ 
nates:  the  equations  no  longer  contain  singularities  at  the  poles,  and  the  primitive  form  of  the 
equations  (momentum  form)  can  be  written  in  conservation  form.  The  conservation  form  of 
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the  equations  are  the  natural  form  of  the  equations  because  they  arise  from  the  conservation 
laws  for  mass,  momentum,  and  energy.  For  this  reason,  this  form  conserves  all  of  these 
quantities  better  than  the  non-conservation  form  and  has  been  observed  to  give  better  results 
than  the  non-conservation  form,  particularly  in  the  presence  of  strong  gradients  [26],  Our 
results  show  that  mass  and  energy  were  conserved  up  to  99.9  per  cent. 

By  mapping  the  three-dimensional  Cartesian  physical  space  to  a  local  computational  surface 
space  we  were  able  to  construct  high-order  spectral-type  finite  elements  on  icosahedral  grids  on 
the  sphere.  The  high-order  basis  functions  used  take  into  account  the  curvature  of  the  element 
as  it  lies  entirely  on  the  sphere;  in  fact,  the  surface  functions  are  the  typical  two-dimensional 
planar  basis  functions.  As  a  result,  although  we  discretized  the  equations  in  three-dimensional 
Cartesian  space,  the  problem  then  reduced  to  the  integration  of  high-order  surface  element 
integrals.  In  addition,  a  new  spectral  icosahedral  grid  has  been  introduced  which  offers  many 
possible  grid  configurations  because  it  is  based  on  the  generalized  icosahedral  grid  introduced 
previously  by  the  author  [12].  The  ideas  presented  in  this  paper  are  generalizable  towards  the 
construction  of  a  parallel  baroclinic  model  (i.e.,  fully  three-dimensional  global  atmospheric 
model),  which  is  the  goal  of  our  current  research. 
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APPENDIX  A.  MAPPING  TO  COMPUTATIONAL  SPACE 

In  this  appendix  we  show  that  the  method  described  in  this  paper  contains  the  Eulerian 
method  described  in  References  [11,12]  for  the  special  case  of  linear  triangular  elements.  To 
summarize  the  method  introduced  in  Reference  [12],  we  give  the  linear  basis  functions  on  a 
triangle  in  three-dimensional  Cartesian  space  as 


1 

W\ 


(a,x  +  bty  +  c,z) 


(19) 


where 


\j\  =Xi-(x2  X  x3) 


(20) 


and 


G  =  y?k  -  y^p  bt  =  xkZj  -  XjZk, 


Ci  =  XjVk  -  XjJj 
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The  indices  i,  j,  k  =  1, .  .  . ,  3  are  the  vertices  of  the  triangle  and  are  cycled  in  the  following 
order: 

k 

/  \ 

i  -»■  j 


and  the  magnitude  of  the  normal  vector  to  the  surface  element  is  given  by 
|G|  =  Kxj-Xj)  x  (x3-x1)| 

In  order  to  show  that  both  methods  are  equivalent  on  linear  triangular  elements,  we  need  only 
show  that  the  matrices  given  by  Equations  (5)— (8)  are  identical  for  both  methods.  Starting 
with  the  simplest  of  these,  let  us  begin  with  Equation  (5).  In  the  case  of  linear  elements,  the 
basis  functions  in  local  space  Q  are 


\-%-ri 

S 

>1 


(21) 


If  we  let  £  =  \j/ 2  and  //  =  iJ/3  then  we  could  write  ^  =  1  —  £  —  tj,  thereby  yielding  the  exact  same 
basis  functions  for  both  methods.  In  addition  by  substituting  Equation  (21)  into  Equation  (13) 
we  get 


G 


dx  dx 

dd  X  dr] 


|(x2-x1)  x  (x3  x3 ) | 


thereby  showing  that  the  mass  matrix  M  is  identical  for  both  methods.  This  also  means  that 
any  matrix  not  having  a  derivative  will  also  be  identical,  such  as  the  rotation  matrices  R'.  Next 
we  show  that  the  derivatives  are  equivalent  in  both  methods. 

From  Equation  (9)  and  the  chain  rule,  we  can  write  the  x  derivative  (where  the  y  and  2 
derivatives  would  result  in  similar  forms)  of  the  basis  functions  as 

di//j  _  8\j/j  d£  d\jjj  dri  r  dtj 

dx  ^  di;  dx  ^  dtj  dx  '  dx 


which,  after  substituting  Equation  (21)  gives 
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<¥,■ 

8x 

The  metric  terms  d^/dx  from  Equation  (12)  are 

dd,  1 

j-x  =  ^  [(>’3  -  y  i)z  -  y(?  3  -  -i)] 

dr]  1 

=  |7|  ^  (“2  ~  ~  ('V2  “ 

dC  1 

^  =  |j|  [( v2  -  v i  )(23  -  zL)  -  (j3  -  jOCzj  -  zO] 

Substituting  these  metric  terms  gives 

=  jjj  [  -  (J3  -  y l)z  +  T(>3  -  -l)  -  T(z2  -  -l)  +  O2  -  Vl)-  +  ^l(j2  -  Tl)(z3  -  -1) 

-  'AlO'3  -  Tl)02  - -1)] 

where  £  =  1  has  been  used.  After  cancelling  terms  we  get 
3ij/  1 

=  jjj  [( V2  -  y3)z  +  y(z3  -  z2 )  +  ^(>2  -  Ji)(z3  -  Zj)  -  ^(>3  -  yx)(z2  -  Zj)] 

Recall  that  we  can  expand  the  co-ordinates  in  terms  of  the  basis  function  approximation 

x=  X  £&x« 

1=  1 

Including  this  approximation  for  y  and  z  above  and  cancelling  like  terms  yields 
c)\J/  1 

-^7  =  pj  [(^1  +  ^2  +  &)(  T2-3  -  y  3^2)] 

However,  the  sum  of  the  basis  functions  (since  they  are  cardinal  functions)  are  equal  to  1, 
thereby  giving 


—  £  t —  Cx~  +  i//j 
OX  ox 

r8£  r  dC 
cTx  +  ^Yx 
d) ]  r  dC 
V.v  'Aw 


dC_ 

dx 
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#1  1  ,  , 
to- R 

From  Equation  (19)  it  is  immediately  obvious  that 

#i  1 

and  so  we  must  now  show  that  |,/|  =  |Jj.  (Note  that  the  same  procedure  for  the  rest  of  the  basis 
functions  results  in  similar  forms.)  From  Equation  (11)  we  can  write  the  Jacobian  in  Q  as 


which,  after  substituting  the  approximation  for  x  yields 
J=  +  lj/2X2  +  ^3X3) -[(X2  -  X,)  X  (X3  -  Xj)] 

and  expanding  the  vector  product  terms  gives 

J  =  (t^  +  lj/2X2  +  ^3X3)  -[x2  x  x3  -  x2  x  x,  -  x,  x  x3] 

Next,  expanding  the  scalar  product  gives 

J=  ^!Xj-(x2  X  X3)  -  lA2X2-(x,  X  X3)  -  \j} 3X3 -(x2  x  xt) 
which  can  be  rewritten  as 

J  =  (A,  +  <j/2  +  ifijx  1  •  (x2  x  x3) 

Since  the  sum  of  the  basis  functions  equals  1,  then  we  get  an  identical  expression  to  Equation 
(20).  Therefore,  any  matrix  containing  a  derivative  is  identical  for  both  methods,  which 
includes  the  remaining  two  matrices  A  and  /JX;  therefore  the  Eulerian  method  given  in 
Reference  [12]  is  contained  in  the  method  described  here. 
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